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THE NUMERICAL CALCULATION OF LAMINAR 


BOUNDARY- LAYER SEPARATION* 

John M. Klineberg and Joseph L. Steger 
Ames Research Center 


SUMMARY 


Iterative finite-difference techniques are developed for integrating the 
boundary- layer equations, without approximation, through a region of reversed 
flow. The numerical procedures are used to calculate incompressible laminar 
separated flows and to investigate the conditions for regular behavior at the 
point of separation. Regular flows are shown to be characterized by an inte- 
grable saddle-type singularity that makes it difficult to obtain numerical 
solutions which pass continuously into the separated region. The singularity 
is removed and continuous solutions ensured by specifying the wall shear 
distribution and computing the pressure gradient as part of the solution. 
Calculated results are presented for several separated flows and the ccuracy 
of the method is verified. A computer program listing and complete ..olution 
case are included. 


INTRODUCTION 


During the past decade, various approximate methods have been developed 
to calculate separated flows by using the boundary-layer equations. The most 
popular schemes have been integral, or moment, methods based on the early work 
of Abbott, Holt, and Nielsen (refs. 1-3) or Lees and Reeves (refs. 4-8). In 
the integral approach, the boundary- layer equations are multiplied by a power 
of u and converted into a system of ordinary differential equations by inte- 
grating across the viscous layer. Regions of attached and separated flow are 
treated similarly because the average convection in the boundary layer is 
always in the downstream direction. 

A second type of approximate method, first proposed by Reyhner and 
Flugge-Lotz (ref. 9) uses finite-difference techniques (refs. 10 and 11). 

This approach uses a forward-marching procedure, with all convective deriva- 
tives set to zero in regions of reversed flow for numerical stability. The 
conservation of momentum and energy is therefore violated in the portion of 
the separated flow bounded by the zero-velocity line, although the errors 
introduced by this approximation are not expected to be significant for small 
laminar separation bubbles. Both the finite-difference and integral methods 
have produced good agreement with experimental data, particularly for 
compression-comer flows and shock-wave/boundary-layer interactions (see the 
review in ref. 1 2). n 

*Presented at the AIAA 12th Aerospace Sciences Meeting, Washington, D. C. 
Jan. 30-Feb . 1, 1974. 



The first finite-difference integration of the complete boundary- layer 
equations through a region of reversed flow was performed by Catherall and 
Mangier (ref. 13) . This report provides the best previous numerical evidence 
of flows that are regular at separation. A continuous solution was obtained 
by specifying the displacement thickness downstream of an appropriate point 
near separation and determining the pressure gradient by streamwise integra- 
tion. The numerical procedure developed instabilities in the reversed-flow 
region, however, and the integration was continued only by decreasing the con- 
vergence criterion at each station. As the authors point out, this difficulty 
is to be expected because the region of separated flow should actually be 
integrated in the upstream direction, with boundary conditions provided from 
downstream. 

There have also been several numerical studies of nonlinear parabolic 
equations of mixed type, where the direction of increasing "time" reverses in 
some region of the flow field. One of these investigations, by Klemp and 
Acrivos (ref. 14), considers the flow over a finite, stationary flat plate 
whose surface moves at a constant velocity opposite that of the free stream 
(i.e., a rotating belt). The boundary layer is divided into two regions along 
the unknown zero-velocity line and the equations are integrated in the appro- 
priate flow directions, with the final solution obtained by iterating for the 
location of the common boundary. It is not evident that this technique would 
prove effective for calculating boundary- layer separation because the region 
of reversed flow results only from the upstream motion of the surface of the 
plate. Also, the pressure gradient is assumed to be zero and the shear 
stresses remain positive throughout the flow field. The singularities at 
separation and reattachment are therefore caused by discontinuities in the 
boundary conditions and are not associated with the vanishing of the surface 
skin friction. 

A more useful numerical procedure for calculating the flow past an 
impulsively started flat plate has recently been developed by Dennis (ref. 15) . 
For this problem, the motion at short times is given by Rayleigh's error 
function solution, while the final steady-state condition is given by the 
Blasius profile. Although the transition from the initial to the final state 
can be calculated directly in the three independent variables (ref. 16), 

Dennis formulated the problem in similarity coordinates where the governing 
equation is parabolic and of mixed type. The convective derivatives were 
approximated by backward or forward differences where appropriate, and the 
solution was obtained through a successive overrelaxation procedure. This 
numerical technique with certain modifications can also be applied to the 
equations that describe boundary-layer separation. The two problems are, of 
course, different in many important respects. In particular, there is nothing 
corresponding to reattachment for the impulsively started flat plate, and the 
downstream (large time) boundary conditions are given. One of the more 
interesting features of boundary- layer separation is that although there is an 
embedded region of reversed flow and of upstream influence, the overall 
problem remains parabolic in the downstream direction. 

The present investigation develops a numerical procedure for integrating 
the laminar, incompressible boundary-layer equations, without approximation, 
through a region of reversed flow. Under Development of Numerical Method, a 


2 



model problem is examined to determine convergence and stability criteria, and 
iterative finite-difference schemes are developed to solve the nonlinear 
equations. Under Results and Discussion, the numerical procedures are used to 
investigate the conditions for regular behavior at the point of separation. 

The separation (and reattachment) points are shown to be saddle-type singu- 
larities in the physical plane, which make it difficult to obtain numerical 
solutions that pass continuously from the attached region to the separated 
region. The singularities are effectively removed, however, by specifying the 
wall shear distribution and determining the pressure as part of the solution. 
These inverse calculations are used to infer the type of pressure distribution 
required for the boundary layer to pass smoothly into a region of reversed 
flow. Where possible, results are compared to relevant analytical or similar- 
ity solutions to verify the accuracy of the calculations. The extension of 
the method to compressible flows and to the solution of complete viscous- 
inviscid interactions is indicated in a separate section. 


DEVELOPMENT OF NUMERICAL METHOD 


The Differential Equations 

The boundary- layer equations for two-dimensional, laminar, incompressible 
flow are 


|H + |X = o 

3x 3y 


3u 3u ^ U e 3 2 u 

u — + v — = u j — + 

3x 3y e dx 9y2 


(la) 

(lb) 


where the Reynolds number has been ejgjlicitly removed by introducing the 
usual scaling x = x*/£, y = (y*/£)A, u = u*/u ro , v = (v*/u m ) SR, R = u ro *£/v. 
Here superscript (*) indicates the physical or untransformed variable. 

Boundary conditions are u = v = 0 and u -* u e as y In a direct 

problem, u e is specified as a function of x, while, in an inverse problem, 
an alternate condition such as (3u/3y)o or v e is given as a function of x. 

In this case, u e must be determined as part of the solution process. 

The parabolic nature of the equations is evident in von Mises coordinates: 


3z 

3x 


3 2 z 

dip 2 


(2a) 


u 


dip 

dy 


(2b) 


with u e 2 - u 2 = z and v = -d\p/dx. Equation (2a) is clearly a heat equation 
in which the coefficient u changes sign in regions of reversed flow. Because 
there is no downstream boundary condition, the solution is not unique unless 
the separated zone is entirely confined within the domain of integration. 
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The equations can also be written as a system of nonlinear first-order 
equations in conservative form, for example. 


_ 9 _ 

8x 



-+) 


3u 

8y 


0 

3y 

(3a) 

o 

9y 

(3b) 

■ g = 0 

(3c) 


Because the equations are nonlinear, discontinuities may occur in the flow 
field even though continuous boundary conditions are specified. Equations 
(3a), (3b), and (3c) possess the following weak solutions: 


[u 2 - u^sin 0 = [v 2 - Vi]cos 0 


(4a) 


[u 2 2 - Ui 2 ]sin 0 = [g 2 - gi ] cos 0 


(4b) 


0 = [u 2 - ujjcos 0 (4c) 

where u e is assumed to be continuous and 0 is the angle between the axis 
and a plane separating conditions 1 and 2. If 6 < tt/2, equation (4c) ensures 
that u 2 = uj and, consequently, all the variables are continuous. When 
0 = tt/2, the weak solutions are indeterminate. In particular, v may be dis- 
continuous with a jump of indeterminate strength even with u continuous. 
Furthermore, if u is discontinuous, then, from equation (4a), [v 2 - vj] 


Preliminary Numerical Considerations 

As equation (2a) in particular shows, in the separated- flow region, 
information must be allowed to propagate upstream with the reversed flow 
velocity. A natural way to fulfill this requirement consistent with restric- 
tions of numerical stability is to treat the x-derivatives with backward 
(upwind) finite-difference formulas in attached flow regions and with forward 
(downwind) finite-difference formulas in the reversed flow region. However, 
this means that at least a portion of the difference equations will require 
simultaneous solution. Furthermore, the extent of the separated region is 
unknown and, because the equations are nonlinear, an iterative finite-difference 
method appears to be the most efficient way to find a solution. Here, of 
course, one can rely on experience obtained with type-dependent relaxation 
methods employed for transonic flow fields (refs. 17 and 18). 

As an alternative to a type-dependent differencing scheme, interpolative 
(elliptic) finite-difference formulas such as central differencing can be used 
over the entire flow region. In fact, in the absence of discontinuities. 
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parabolic and hyperbolic problems can be solved with interpolative 
differencing, provided the boundary conditions are properly satisfied. Of 
course, for a simple initial -value problem, a marching process that uses back- 
ward differencing is generally far more efficient than a simultaneous solution 
process. 

The choice of whether to use backward -forward differencing, central 
differencing, or some hybrid of these will depend on the efficiency and 
accuracy obtainable in the iterative finite-difference method. In any case, 
no downstream boundary conditions can be supplied for the boundary- layer 
equations, so the last computed profile must be attached to allow the use of 
backward differencing for the x-derivatives . 

The success of a numerical method also depends on the choice of 
variables into which the equations are cast. Equation (2), for example, is 
not suitable because the variable \p is multivalued in the separated region. 
For the most part, equations (la) and (lb) appear to be the most appropriate 
to difference with a high probability of being readily extended to more 
complex (e.g., three-dimensional) flows. 

Because the boundary- layer exhibits extensive growth in the. x-direction, 
it is essential for numerical efficiency that this growth be scaled out. 

This can be accomplished by introducing a variable, growing grid system or 
by using a transformation that keeps the viscous layer of nearly uniform 
thickness. The following transformation is used: 



so that equations (la) and (lb) become 

3u 3v m + 1 - n 

x — + — + — = — u = 0 (6a) 

3x 2 

3y 

xu |^- + v — = m(l - u 2 ) + — — (6b) 

3x 3y »y 2 

Boundary conditions are indicated in figure 1. These equations can also be 
written as a single equation for the stream function 

f"» + ff" + m(l - f' 2 ) = x(f’f ’ - ff") (7) 

4 A X 
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Figure 1.- Difference operators and boundary conditions for 

relaxation calculation. 


Iterative Finte-Difference Method 

In the first stages of developing a finite-difference method, it is 
useful to begin with the study of a model problem. A model equation is 
obtained here by linearizing equations (la) and (lb); the iterative convergence 
criteria are reviewed and an appropriate choice of difference formulas is made 
so that the simple model equation is iteratively stable. In the following 
section, the convergence of the difference equations to the differential equa- 
tion is considered, and iteratively convergent differencing schemes for the 
nonlinear boundary- layer equations are subsequently given without analysis. 


Model problem- Equations (la) and (lb) are simplified with 

u = Uq + u 
v = vo + v 

so that the model equation becomes 



( 8 ) 


3 2 u 3u 3u 

9y2 - u o 9x - v 0 ay 


b (x) 


(9) 


In any local domain, uq and v 0 are treated as constants. Equation (9) also 
represents the transformed equations, equations (6a) and (6b), if an average 
value for xuq is substituted for uq. 


If convergent difference algorithms and convergent iterative procedures 
can be selected for equation (9), subject to all reasonable choices of u 0 
and vq, it is assumed that such schemes can be successfully adapted to equa- 
tions (la) and (lb). While explicit and implicit marching procedures have 
been developed and extensively studied for parabolic equations of standard 
type, a comparable development does not exist for relaxation schemes. The 
development of such a scheme is undertaken below where the primary concern is 
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to ensure that the relaxation procedure is valid for both positive and 
negative values of uq. 

Iterative convergence criteria - Once equation (9) is differenced over a 
discrete network of grid points, one is left with the task of inverting the 
linear system of equations 


Au - c = 0 


( 10 ) 


where the components of u consist of the dependent variables at each grid 
point. Then the most general first-degree iteration scheme for equation (10) 
is 


. ;("> = hH[A :' n ) - t] 


0D 


where H is a conditioning matrix usually implicitly built into the iterative 
solution algorithm; here we chose to extract a parameter h from H. It 
should be understood that this type of iterative solution algorithm can treat 
nonlinear equations with the same ease as linear equations. 

Equation (11) has the recursive solution: 



(I + hHA) n u ( -°' ) 


n-1 

^ ^ (I + hHA) m hHc 
m=0 


( 12 ) 


so if the matrix (I + hHA) has a spectral radius (i.e., largest eigenvalue in 
absolute magnitude) less than 1, then (I + hHA) n -*• 0 for n sufficiently 
large. Furthermore, from the Neumann lemma (ref. 19, p. 26, or ref. 20, 
p. 82), it is evident that 


n- l 

-£ 


m=0 


(I + hHA)^ 


- A" 1 


(13) 


or u -*■ A -1 c as required. Thus the sufficient condition for iterative 
convergence is that all 


V E 


1 + 


ha . I 
J 


< 1 


(14) 


HA. Hence, if all the real parts of the 


where Oj are the eigenvalues of 

c r i Vil r* Amrs 1 qv rf • a-P V* « 

j jxl/ij v- x -j ux v j_ cnv OCilll^- 

vergence. This is an asymptotic convergence criterion for n sufficiently 
large. For the scheme defined by equation (11) to be efficient, the matrix 
HA must not have a large condition number (refs. 19 and 21) nor should the 
imaginary parts of aj be large compared to the real parts. The eigenvalue- 
convergence criterion does not guarantee that the norm of (I + hHA) n u(°) will 
not grow appreciably during intermediate iterations - a situation likely to 
occur if the eigenvectors of HA are linearly dependent or almost so. 
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Convergence of the model problem- The advantage of studying the model 
problem is that analytic expressions are obtained to describe its behavior for 
various choices of differencing. It is assumed that the nonlinear problem 
will share at least some common features. Here, let 


3 z u 




jk 




jk + 


0 


+ 0(Ay 2 ) 


(15a) 



(15b) 


and 


3u 

3x 


jk 


■ ( i )^ 1 




2au 

jk 


(1 - 




+ 0(Ax) (16) 


where a = 1 is first-order backward, a = 0 is second-order central, and 
a = -1 is first-order forward. Using these approximations in equation (9) 
with 3 = Ay/2 and y = (Ay) 2 /2Ax, one obtains 


(1 + v 0 e)u jk _ 1 - 2u. r + (1 - v 0 3)u. k+i + (1 + a)yu 0 u ;j _ ijk - 2ayu 0 u. k 

- (1 - a)yu 0 u j+1>k = b. 


(j = 2, 3, 4 


max 


2, 3, 4 . 


K - 1) 
max 


(17) 


If u is the vector whose components are the Ujj. over the ordered grid 
points, equation (17) can be written as the linear system of equations, 
equation (10). The eigenvalues of' A are given by 


(k = 1, 2 . 


a jk = _2 [ 1 + A 1 + v 0 8)(l - v 0 3) cos ( K - +- 1 ) ] 
-2u 0 yja + / -(1 + a) (1 - a) cos (jTt)] 


, K; j = 1, 2 . . ., J; K = K - 2; J = J - 1) 
’ 5 J ’ ’ ’ max max 


( 18 ) 
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where u is assumed to be given on a boundary as needed. If a is 0 or 1 
when up > 0 or if a is 0 or -1 when uq < 0, the cr roots always have 
negative real parts and A is a stable matrix. Thus the point-iteration 
scheme with H = I and h = to/ (2 + 2 uqyo 0 is proven to be convergent for an 
appropriate to < 1 . As another example, the point successive overrelaxation 
(SOR) method has the roots 

^1 - -j^(l + ayuo) = A - oJ -,1 + v 0 S)Cl - v 0 3) cos (^T~t) 

-u 0 y / (1 + a) (-1 + a) c ° s (/f T 

(j = 1, 2 . . . K; j = 1, 2 . . . J) (19) 

and is also iteratively convergent with h = -1 and a proper choice of the 
relaxation parameter to. 



Equations (18) and (19) show that the roots will be complex if a = 0 or 
if | v o 3 1 > 1. This can be detrimental to the convergence rate of a first- 
degree iteration scheme if the imaginary parts become large enough, so the 
central differencing should be restricted to regions where uqy is small. 

The product Vq 8 is normally expected to be less than 1 in absolute value and 
thus has the beneficial effect of reducing the term cos kir/(K + 1). 

In place of the complex roots that occur for a = 0, when a or -1, 
the eigenvectors of HA appear in multiples of the number of J grid points. 
Under these conditions, the norm of an iteration matrix can be expected to 
grow before it decays; however, study of anc * induced matrix norms 
(ref. 19) for the point iteration scheme shows that residual growth cannot 
occur if the spectral radius is kept less than 1. Conversely, numerical 
experimentation with the heat equation demonstrates that the SOR forward- 
differenced scheme (a = -1) swept from left to right can experience appreciable 
residual growth if Ax << (Ay) 2 . If swept from right to left, the residuals 
decay rapidly. 


Convergence to the Differential Equations 


Although the previous analysis shows that iteration algorithms can be 
used to find a solution to the system of difference equations, it does not 
prove that the solution of the difference equations will converge to the solu- 
tion of the differential equations as the grid is refined. However, with the 
exception of the central differencing scheme, all the schemes to be introduced 
are known to be stable and consistent for the heat equation (cf . ref. 22) . 


If one assumes periodic boundary conditions in x and end conditions in 
y, then sufficient conditions for convergence of the centrally differenced 
heat equation 
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are Ax - 0(Ay 2 ) and Ay - 0(Ax 2 ). (This is not an explicit leap-frog scheme.) 
Here convergence implies that the difference between the exact solution to 
the differential equation and the exact solution to the difference equation 
will vanish as the grid is uniformly refined over a fixed domain. That is, 
the summation of truncation errors given by A _1 e 0 as Ax, Ay 0 where_^ A 
is the matrix formed by the difference equations over both y and x, and e 
is the vector of truncation errors. While the complete convergence proof is 
too lengthy for this report, note that A is a normal matrix and hence is 
unitary similar to a diagonal matrix of its eigenvalues (ref. 21) . The 
eigenvalues are 

c. k ( 2 4xA) = 4 [l - cos (j^+"t)] + 2H±Dsi"(nr) 


(j = 1, 2 . . . J; k = 1, 2 . . ., K) (20) 

and || A -1 1| 2 = (min |tfjjJ) -1 s0 1 1 A” 1 [ | 2 Il e ll2 i s simply determined. 


Finite-Difference Equations and Solution 


Two second-order-accurate differencing schemes were developed to solve 
the boundary- layer equations (6a) and (6b) . The first of these proved 
superior for the separated flows computed in this investigation. The second 
more conventional method is described because it may prove efficient for 
certain extensions of the present approach. 


The first method employs the central-differencing schemes for Uyy and Uy 
given by equations (15a) and (15b) . The term xuu x in equation (6b) is 
. regrouped as 0.5x(u 2 ) x and backward-differenced: 


x 8u 2 
2 9x 


jk 

(B) 




- 4u? , + 

3-1, k 


2Ax 



+ 


0 (Ax 2 ) 


( 21 ) 


for u > 0.015 or j = J max . When u < 0.005 or if j = 2, central 
differencing is used: 


x 

2 


Ini 

3x 


jk 

(C) 



0 (Ax 2 ) 


( 22 ) 


10 



In the intermediate zone, 0.005 * u * 0.015, the backward and central formulas 


3uf 

3x 






n 

1_ 

jk 2 

,, ,, 3u 2 

(1 * 3x 

+ (1 
jk 
(B) 

- a) 9u2 
} 3x 

n 5^ 

v — / 

1 


(23) 


with a = 1 + 200(u - 0.015). The difference stencils are indicated in 
figure 1. 


We emphasize that the blending defined by equation (23) is used solely to 
enhance the iteration process and is not otherwise fundamental. It is obvious 
that when the difference equations are switched at a given value of u', a dif- 
ferent set of data points is sampled and slightly different truncation errors 
result. The change in the residual error vectors at this point can be large 
enough to drive u( n+ l) back across the value at which switching occurs. 

This can then start an oscillatory mode with little decay. The blending sim- 
ply modifies the differencing relations in a continuous fashion so that the 
residuals vary smoothly. In the present scheme, the blending is completed at 
u = 0.005 to avoid a special operation at separation and reattachment. The 
blending can also be used between 0 < u < 0.01 without changing the results. 


The continuity equation is differenced with the modified Euler scheme 
(i.e., trapezoidal rule or Crank-Nicholson differencing): 


v,,. - v 


jk jk-1 2 y3y 


/ 3v 


3v 

y3y 

jk 

+ 3y 


with 


9v 



L 

2, J - 1 

max 

. , and 

3v 



3y 

jk 

x 

LA 

f max' 

Note 

that 


jk-i) 

( u. i -u. , ,\ m. +1 

lzlA) + _J u 

2Ax J 2 jk 


+ 0 (Ay 2 ) 


+ 0 (Ax 2 ) 


ji 


- 4u . , , + u . „ , \ m. +1 
— izlA izlA) + _2 u 

2 Ax J 2 u jk 


(24) 


(25a) 


+ 0 (Ax 2 ) (25b) 


at j — HUIC 1.110.1- AU^ 

J max ) in both the attached and reversed flow regions. While equation (24) is 
generally recommended, two schemes implicit in the y direction are presented 
as alternatives. Either the second-order-accurate "shifted" scheme 


^ + 0(Ay 2 ) (26) 

(where point jk is updated in the relaxation) or the third-order accurate- 
in-y "abated Hermite" scheme: 
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3v 

3y 


jk-l 


£ 

+ 6 3y 


jk 


£ 

6 3y 


jk+l ) 


+ 0(Ay 3 ) (27) 


can be used with Sv/Syh^ again defined by equation (25). Both alternative 
schemes generate diagonally dominant tridiagonal blocks if a backward two- 
point differencing is used at the edge where v varies linearly. Effectively, 
equations (24), (26), and (27) give the same results. 

An additional difference algorithm must be introduced if an inverse 
problem is solved. To impose a specified shear distribution, the momentum 
equation is evaluated at the surface: 


m 


3 2 u 



(28) 


The second derivative is differenced as a function of x(x) to generate the 
second-order-accurate relation: 


3 z u 


3y' 


jl 



+ 8u. „ 

J2 


2 (Ay) 2 



fill 

Ay 


(29) 


Wake flow is treated in the same fashion with x = 0 and the centerline 
velocity Uq specified: 


(1 - u 0 2 )m = - 



du 0 

+ XU ° dT 


(30) 


With the choice of differencing established, the solution procedure is 
straightforward. An approximate solution is input, usually by assuming a 
Blasius profile with m = 0 everywhere. For an inverse problem, a new distri- 
bution of m is then predicted for the specified boundary condition using 
either equation (28) or (30), with m updated by the relaxation (here 
written for x(x) specified): 


m. 

J 


(n+1) 





8 u i 2 + u p \ 

2 (Ay) 2 / 


(31) 


For a poor guess of the initial solution, to is initially kept small, 
a) = 0(0.05). New values of u are then found from relaxation of the momentum 
equation, while new estimates of v follow from continuity. This iteration 
sequence continues (with u increased as the initial guess is improved) until 
an equilibrium or converged state is reached. 

Solutions are found by both point and line successive underrelaxation 
(SUR) by using the iterative correspondence: 
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+ u)H R 
u u 


' (32a) 


-Kn+l) ^(n) 

u l J ' 


-»-(n+i) 
v 1 



+ u)H R 
v v 


(32b) 


The residual vectors R u and R v represent the differenced momentum and 
differenced continuity equations, H u and H v are the conditioning matrices of 
the SUR algorithm, and uj is the relaxation factor. The line method (not 
used in eq. (32b)), in general, converges faster than the point scheme, but it 
is more sensitive to changes in m, making it more difficult to control in a 
computer batch mode. For moderate reversed flows and grid spacings with Ax 
approximately equal to Ay, the optimum relaxation parameter is slightly less 
than 1 for point SUR with co = 0(0.5) for equation (31). For line SUR, the 
optimum relaxation parameter is 0(0.4) and a> is the 0(0.15). The point SUR 
method fully converges in 400 to 800 iterations for a grid of 80 j -points and 
50 k-points. Highly separated cases with rapid variations in the flow 
quantities require the higher iteration counts. 

Note that, when u is negative, it is possible to blend from the 
central into a three-point forward difference and that this variant of the 
relaxation procedure is iteratively convergent. For very large reversed-flow 
regions, it may be advantageous to program this additional logic. Experience 
also shows that switching at u = 0 from a three-point backward differencing 
into a three-point forward differencing without blending first into the 
central differencing is not convergent. 


The second method developed is patterned after the Crank-Nicholson scheme. 
Equations (6a) and (6b) are first put into conservative form 


3v 3xu (m - l)u _ 

7T + a3T + -Hr — = 0 
3y 


(33a) 


3 fuv) 

ay 


a 2 u 3(xu 2 ) 3m - 1 ~2 
+ _ — l. + 

3y 2 3x 2 


m = 0 


(33b) 


The continuity equation is treated as before, and the y derivatives in the 
momentum equation are again centrally differenced by use of relations (15a) 
and (15b) . The x-differencing is Crank-Nicholson 


(xu 2 ),. - (XU 2 ), 

JK J ~ 


9 (xu 2 ) 


jk 


3 xu 2 

3y 


Ij-i.k 


(G > 0.01 or j = J max ) 


with 


3 (xu 2 ) 


3^u 


3x 


jk 3y 2 


+ m. - 


jk 


i - hr-h 


3uv | 
jk " 3y 
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(34) 


(35) 
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where the appropriate central-difference formulas are substituted for the y 
derivatives. For reversed flow, forward differencing is used 


(xu 2 ) 


j + lk 



Ax 

8 (xu 2 ) 

2 

9x 



(u < -0.01) (36) 

and the two schemes are linearly blended in the interval -0.01 < u < 0.01 (see 
fig. 1). As before, the blending is used solely to enhance the iteration 
process. 

The Crank-Nicholson scheme has been solved by both point and line SUR, and 
for either process the relaxation parameters are approximately those described 
for the previous point method. This second method requires slightly more 
algebra per step and, in general, has a slower rate of convergence than the 
first method. 


The conservation- law form of the Crank-Nicholson method is not 
considered to be an advantage, and the procedure generally predicts m dis- 
tributions that are slightly oscillatory. The oscillations decay as Ax/ Ay 
decreases, and they are. confined to the relatively uninteresting attachedflow 
regions. Of course, m is a sensitive function of the solution and the u and 
v distributions are much smoother. A nonconservative version of the Crank- 
Nicholson scheme was also programmed. In this case, the oscillations in m 
were negligible in attached-flow regions but observable in the separated zone. 


Finally, we remark that a very stable first-order-accurate method can be 
developed by replacing the x differencing by 


and 



(37a) 


(37b) 
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This scheme, with xu x of continuity also first-order-accurate and switched in 
an identical fashion, will generally give "computational results" for the 
first problem, m specified. This first-order method is not recommended 
because a much finer x-grid spacing is required to maintain accuracy. This 
method proved useful for the numerical experiments described in the next 
section. 


RESULTS AND DISCUSSION 


In this section, the iterative finite-difference procedure is used to 
integrate the boundary- layer equations through a region of reversed flow. The 
separation -point singularity is investigated and conditions for regular 
behavior are determined. Calculated results are presented for a number of 
separated flows and the accuracy of the method is verified. Possible indica- 
tions of the breakdown of the boundary -layer assumptions are also examined. 


Direct Solutions 

One of the most extensively studied problems in separating boundary- layer 
flows is the response of a flat-plate boundary layer to a linearly retarded 
external stream. This problem has been investigated by Howarth (ref. 23), 
Hartree (ref. 24), and many others; recent solutions have also been obtained 
by Briley (ref. 25) and Leal (ref. 26) using the full Navier-Stokes equations. 


A sequence of calculations for different mesh spacings is shown in 
figure 2. The external velocity was specified to decrease linearly from the 



X 


Figure 2.- Calculation for linearly 
retarded flow. D 


origin and the first-order-accurate 
difference scheme was used because 
of stability considerations. As the 
mesh is refined, the separation 
point moves upstream, with the last 
calculation in exact agreement with 
the accurate results of Smith and 
Clutter (ref. 27) for the flow 
upstream of separation. In the limit 
of zero-mesh spacing, it is evident 
that the solution is singular, with 
the wall shear approaching zero as 
the square root of the distance 
upstream of the separation point and 

Lilt? I lO 1'llttl 1 VtiUtiLy uuuuuaucu. iii-Lo 

type of behavior has been discussed 
in detail by Goldstein (ref. 28) and 
by Landau and Lifshitz (ref. 29) 
among others. 


The interesting result here is that the use of an iterative finite- 
difference scheme which contains type-dependent operators allows the solution 
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to be "continued" in the downstream direction. As the mesh is refined, it 
becomes evident that the flow fields upstream and downstream of separation are 
essentially independent, and the solution is therefore not meaningful. The 
wall shear jumps discontinuously to a negative value at separation and the 
normal velocity v becomes unbounded; all flow quantities subsequently remain 
continuous downstream of the jump and through reattachment. The magnitude of 
the discontinuity is determined by the specified pressure distribution in the 
separated zone. In a set of simple numerical experiments, a constant external 
velocity distribution was smoothly joined to the linearly retarded flow at 
different values of x. As the joining point was moved downstream, .the magni- 
tude of the jump and the extent of the reversed-flow region increased 
monotonically, with separation remaining at x = 0.96. 

Singular behavior at the point of separation is thus related to the fact 
that the wall shear x = (9u/3y)o is nonanalytic; in particular, 
x ~ (x s - x) V 2 and 9x/9x ■> ® as x + x s , the separation point. Therefore, the 
most obvious means of ensuring regular solutions at separation is to specify a 
continuous wall-shear distribution. The pressure distribution can then be 
determined as part of the solution by satisfying the momentum equation at the 
surface. Note that because the equations are nonlinear, it is not possible to 
guarantee that discontinuities will not occur in the flow field even with 
analytic boundary conditions prescribed (see ref. 30 for hyperbolic equations, 
or the weak solutions, eqs. (4a), (4b), and (4c)). 


Inverse Solutions 

With the wall-shear distribution specified, m can be determined from 
equation (31) and the second-order-accurate differencing scheme generates con- 
tinuous solutions that give no indication of singular behavior at either 
separation or reattachment. These solutions are demonstrated to be regular 
under Accuracy Check. An inverse calculation cannot be duplicated by the 
direct method, however. Starting with a fully converged inverse solution, the 
calculation diverges if the iteration is continued with m fixed, that is, 
the relaxation parameter u is set to zero in equation (31) . Two examples of 
this type of inverse (x specified) and direct (m given) calculation sequences 
are shown in figure 3. After as many as 500 iterations (less if the solution 
is initially perturbed), the residuals begin to grow and the relaxation pro- 
cedure either becomes unstable or converges to a different "solution" of the 
difference equations. As the mesh is refined, the second-order scheme fails 
to converge while the first-order method, for moderate grid spacing, generates 
computational results containing a discontinuity. 

The fact that the direct calculation fails to duplicate a converged 
inverse solution -cannot be ascribed to instabilities in the numerical scheme. 
The only difference between the two calculations is the value of the relaxa- 
tion parameter <o in equation (31), and the solution processes are essentially 
identical. The numerical evidence therefore strongly suggests the existence 
of a saddle-type singularity at the separation point. Because of this critical 
point, roundoff and residual errors are sufficient to cause a completely 
converged solution to diverge when the pressure-gradient parameter is held 
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X 

Figure 3.- Inverse/direct calculations 
that indicate existence of saddle 
point. 


where 


fixed. There are no other possible 
sources of error in the calculations: 
the variation of m is determined to 
arbitrarily high accuracy by the 
inverse solution, and no interpolation 
or differentiation is required as for 
computations with experimentally 
determined pressure distributions. 

With the pressure gradient correspond- 
ing to a completely regular flow field 
prescribed, the equations contain a 
saddle-type singularity at separation 
that makes a continuous numerical 
solution difficult to obtain. The 
saddle point is removed from the 
domain of integration, however, by 
specifying the wall shear rather than 
the pressure gradient as a boundary 
condition. A discussion of the 
essential differences between the two 
types of calculations is presented 
below. In the following section, the 
conditions for regular behavior at 
the point of separation are examined. 

Saddle Point 

The difference between the 
direct- and inverse-calculation pro- 
cedures can best be illustrated by 
examining the boundary- layer equations 
near the surface. Expanding the 
velocity profile in a Taylor series 
in y yields 

(38) 


2 3 

u (x , y) = fjy +.f 2 fr + f 3 V* + 


f 3 = 0 and the notation is used 


*1 


T, 



Either fi or f 2 (but not both) is prescribed and all other 
determined as functions of x by the differential equations, 
must satisfy the following set of relations: 


are 
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f 4 - fjf { = 0 ' 

f 5 - 2fi£' = 0 

f 6 - 2f 2 f£ = 0 

£ 7 - 4f_if4 + Sf 4 fJ = 0 • 

f 8 - 5£]i£j - 9£ 2 £4 + 5f 4 f* + 9£ 5 f* = 0 


( 39 ) 



Figure 4.- Series expansion for 
similar separation profile. 


where the prime denotes differentia- 
tion with respect to x. One of the 
fi is given by the outer boundary 
condition that u ■> u e as y + “ (see 
ref. 28). The validity of the expan- 
sion procedure near the separation 
point is demonstrated in figure 4 for 
the particular case of similar flow 
with m = -0.09044, corresponding to 
zero shear (see eq. (7)). For this 
case, the only nonzero coefficients 
multiply terms of the order 
T1 t tn-2 ( n = 2 , . . .), and the 

expansion has been continued through 
the twenty- second power of the normal 
coordinate. 

Direct calculations- For the 
pressure gradient specified, the 
coefficients in equations (39) must 
be determined by integrating the fol- 
lowing system of first-order 
differential equations: 


fifi = *4 

4f! 2 f 'J = fjf 7 + 5f 4 2 

14f 1 3 f 7 = 2f ! 2 [f x 0 + 8 P X ( 2 P XX “ 5 P X P XX )] 
+ 33f 1 f 4 f 7 - 35f 4 2 


(40a) 


(40b) 


/ 

The remaining coefficients are given by the algebraic relations: 
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£ 5 = 2 Pxx £ l 

£ 6 2 ^x^xx 

4f 1 2 f 8 = 9p x (fif 7 + 5f 4 2 ) + 4£! 2 (10p xxx f ! 2 - 13p xx £ 4 ) 
f 9 - 8fl(5p x p xxx - 2p 2 x ) 

If one arbitarily terminates the expansions at this point and assumes that 
flO can be correctly specified, then equations (40) and (41) provide rela- 
tions for all coefficients of the lower-order terms. Given the velocity pro- 
file at a particular station, standard numerical techniques can be used to 
integrate equations (40a) and (40b) to determine the adjacent profile provided 
fl is nonzero. As f^ 0, however, the solution becomes increasingly sensi- 
tive to the calculated value of fi, and numerical errors are propagated in 
the direction of integration. The equations are highly nonlinear, with the 
coefficient fj of the derivatives determined by f 4 , which in turn depends 
on f 7 , f 10 , etc., and on the outer boundary condition. Even for a pressure 
distribution corresponding to a regular solution at separation, the numerical 
integration of equations (40a) and (40b) is unlikely to result in values of 
fl and f 4 that vanish simultaneously. In that event, either f{ will be 
infinite, leading to a square-root singularity, or fi will remain positive 
and the calculation will fail to show boundary- layer separation. 

We emphasize that with the pressure gradient specified, the nonlinear 
equation for the wall shear (eq. (40a)) is inherent to the system of differ- 
tial equations. Even with special procedures that would guarantee that f 4 
vanishes at x = 0, the saddle point would remain to confound the numerical 
solution process. The behavior shown in figure 3 is to be expected because a 
converged solution perturbed by small roundoff and residual errors cannot 
remain converged in the presence of the saddle-point singularity. 

Inverse calculations- For the wall shear specified and the pressure 
distribution determined as part of the solution, a different system of 
ordinary differential equations results: 

2xf£ = f 5 (42a) 

10 xf 5 = 2 f 8 + 23x x f 5 - 18f 2 (xT xx + x x 2 ) 

40xfs = 5fn + 93x x f 8 + 3f 5 (160TT xx - 201 x^) 

+ 27f 2 (19x 3 - 16xx x - 20x 2 x ) 

x x xx xxx 


including the algebraic relations 




19 



fk = TT 
H X 

t £ 6 = f 2 £ 5 
£7 = t(4tt - t 2 ) 

Tfg = 4(f 2 £ 8 - f 5 2 ) + 26T x f 2 £ 5 - 36f 2 2 (TT xx + T x 2 ) 
x 2 fio = 4f 2 (f 2 £ 8 - f 5 2 ) + 26x x £ 2 2 f 5 -36f 2 3 (TT xx + t^) 

+ t 3 (27t 3 - 24tt t + 28t 2 t ) 

x x xx xxx 

These equations are linear, with the coefficient of the derivatives r 
specified as a function of x. The system is therefore less susceptible to 
numerical error, and although the matrix of coefficients still vanishes at 
x = 0, the saddle-point singularity has been effectively removed. If the 
numerical integration is accurate enough to ensure that £ 5=0 when t 
vanishes, the solution will pass smoothly through the separation point. 

The basic difference between the inverse and direct problems is that, 
for the pressure gradient prescribed, the unknown shear distribution is 
determined by a nonlinear equation that contains a saddle-type singularity at 
separation. For the wall shear specified, on the other hand, the pressure 
gradient is given by a linear equation that is much less sensitive to numeri- 
cal error. This is probably also the case when the displacement thickness is 
prescribed (see ref. 13). The fact that most numerical evidence indicates a 
singularity at separation is therefore misleading because of the difficulty in 
numerically integrating through the saddle point. Of course, not all pressure 
distributions admit a regular solution (as discussed in the following section) . 

An interesting point is that, provided the correct numerical procedures 
are used, no difficulties are encountered at reattachment (see fig. 2 or 3). 

The reason for this is that any numerical errors made at the reattachment 
point are either integrated out of the downstream boundary or upstream toward 
separation. The direction of the flow, and therefore the differencing scheme, 
results in a solution process that allows integration away from the saddle 
point at reattachment but that requires integration into the singularity at 
separation. 

Several numerical experiments were performed to verify these 
conclusions. In one set of computations, the velocity profiles at separation 
and in the immediate vicinity of that point were held fixed after converging 
the inverse calculation. For these cases, the inverse and direct procedures 
gave identical solutions. Similar results were obtained when an artificial- 
viscosity term equal to eu xx was introduced into the difference equations. 

As the coefficient e was decreased, however, the direct calculation would 
again diverge from the inverse solution. 



20 



Pressure Gradient at Separation 


As shown in the previous section, the existence of a regular solution 
requires that £4 = 0 at the point of separation (see also refs. 31 and 32) . 

The coefficients f 5 , f-j , and fg must also vanish at the point of zero shear, 
and the pressure gradient must therefore satisfy certain specific conditions 
to permit the flow to pass smoothly through separation. The constraints on 
the pressure distribution cannot be determined directly because of the saddle 
point, but must be obtained from the inverse, or shear-specified, calculations. 

It is reasonable to expect that only certain pressure distributions will 
admit regular solutions. The separation profile, for example, is determined 
by both the upstream and downstream flows so that some compatibility relation 
must be satisfied at this station. Also, from kinematic considerations, the 
boundary- layer approximation to the vorticity transport equation is 


3m ^ „ 3m 
U 3x 3y 


3^m 

3y 2 


(44) 


where m = 3u/3y and p x = 3m/3y at y = 0. The restriction on the pressure 
gradient at separation can thus be interpreted as a constraint on allowable 
boundary conditions: the normal gradient of vorticity at the surface is 

required to satisfy some local condition for the vorticity to remain continuous 
at the singular point. 

From physical considerations, a constraint on the allowable pressure 
gradient implies that the interaction between the inner viscous layer and the 
outer fluid essentially determines the conditions at separation. Prandtl 
(ref. 33) recognized this in 1938 when he stated that the pressure field 
could not be chosen arbitrarily for the flow downstream of separation "to agree 
with observation." Most numerical solutions of the Navier-Stokes equations, 
including the recent investigation by Leal (ref. 26) in particular, also 
indicate that, when the interaction with the outer flow is included, there is 
no evidence of singular behavior at separation. 


Because of the nonlinearity of the boundary- layer equations, it is not 
possible to determine the precise pressure-gradient condition that permits a 
regular solution. Certain restrictions on the pressure distribution can be 
inferred from the Taylor series expansion and from the numerical solutions, 
however. The acceleration of a fluid particle near the surface, for example, 
can be approximated as follows: 


3u 

3x 


+ v 


3u y 2 
•5— = fs- TT 
3y 2 > 


2rp l 
*xx 3 


P P 
r x xx 


~r~ + T(4TT 
6 ■ xx 


2 ) 

' ; 60 




(45) 


Immediately upstream of separation, r and p x are positive and t x is nega- 
tive. As t 0, the fluid in a stream tube near the surface continues to 
decelerate, and the streamlines continue to move away from the wall provided 
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( 46 ) 


P (v ^7 _+2t ^^ < I tt I + |tt 2 | ^ + . . . 

^xxrx 6 3 J 1 x 1 1 x 1 60 

For the flow to separate smoothly, then, a restriction on the pressure field 
is that 


p <0 as t 0 (47) 

xx 

There will therefore be an inflection on the pressure distribution upstream of 
the separation point. This requirement is consistent with experimental evi- 
dence, and the existence of a "knee” in the pressure curve is often taken to 
indicate boundary- layer separation. 

The numerical evidence suggests that this condition is not sufficient, 
however. All regular solutions, in fact, satisfy the requirement: 

^ > o at r = 0 (48) 

This is a more restrictive condition than that given by equation (47) because 
m x can be negative for p xx negative. The linearly retarded flow considered 
under Direct Solutions, for example, satisfies equation (47) but not equation 
(48). In a series of papers, Meksyn (refs. 34 and 35) has contended that the 
existence of a minimum in m x was a necessary condition for regular separa- 
tion. He cited Schubauer's (ref. 36) measurements of the flow over an elliptic 
cylinder as experimental verification of this requirement. Similar arguments 
have also been advanced as a result of the use of approximate methods to cal- 
culate supersonic viscous-inviscid interactions (see, e.g., ref. 37). 

The most useful means of examining the numerical results is in the 
x - m phase space (fig. 5). Several typical computations are presented, 
including the locus of solutions for similar flow. In these coordinates, x 
is a parameter that varies along the curves, with Ax for the similarity 

solutions. For this limiting curve, dm/dx, and therefore dm/dx, is zero at 
the point of zero shear. All nonsimilar trajectories, on the other hand, have 

positive m x at both the separation and reattachment points. This condition 

was never violated in approximately 30 different calculations using various 
specified shear distributions. Note that the locus of similar flows is some- 
times taken to indicate singular behavior at separation because 
x ~ (mo - m) 1 / 2 and dx/dm + ” at x = 0. The similarity solutions are obtained 

for m x = 0, however, and the limiting value of dx/dx (= m x dx/dm) must be 

carefully determined if an actual flow is replaced by a sequence of similar 
flows. In any event, the condition for regular separation, that m x - 0 at 
the point of zero shear, is satisfied by both the similar and the nonsimilar 
flows . 
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Figure 5.- Phase space representation. 


- m r , m t 

0.0826 - 0.0662 0.0783 - 0.0776 

0.0835 - 0.0660 0.0767 - 0.0769 



Figure 6.- Comparison with similax 
solutions . 


Accuracy Check 

The phase-space representation 
of solutions presents an opportunity 
to verify the accuracy of the numeri- 
cal procedure. The points labeled A 
and B in figure 5, for example, have 
the same value of t and m as a cor- 
responding similarity solution. The 
left-hand side of equation (7), which 
is completely determined by 
t = f"(0) and m, is therefore zero. 
The local x variation vanishes and 
the similar and nonsimilar profiles 
must be identical at those points. 

The velocity profiles calculated by 
the present scheme are compared to 
adjacent solutions of the similarity 
equation (obtained by fourth-order 
Runge-Kutta integration in ref. 38) 
in figure 6. There are essentially 
no differences in the results 
obtained by the two methods. 

With a continuous shear 
distribution specified, the solution 
is constrained to be regular at both 
separation and reattachment. This 
result can be verified by comparing 
the calculated streamline pattern 
with the local solution of the 
Navier-Stokes equations obtained by 
Oswatitsch (ref. 39) (see also Dean 
(ref. 40) and Legendre (ref. 41)). 

At the point of zero shear, a regular 
solution of the Navier-Stokes equa- 
tions requires that the angle of the 
dividing streamline be proportional 
to the ratio of the x derivative 
of the shear and the pressure 
gradient. In the transformed vari- 
ables, the precise condition is 



where 0 is the angle of the dividing streamline. For a prescribed shear 
distribution, the calculated values of m can be integrated in x to obtain 
u e . The flow in the vicinity of separation and reattachment for a refined- 
mesh calculation (Ax = Ay = 0.1) is compared with equation (49) in figure 7. 
The calculated results agree exactly with the local Navier-Stokes solution at 
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the point of zero shear, again demonstrating that the boundary- layer solution 
is regular. 



Figure 7.- Detailed flow field in the vicinity of x = 0. 


Flow-Field Solutions 


As previously mentioned, a number of different shear distributions were 
specified in an effort to determine the behavior of the boundary- layer equa- 
tions in separated flow. Some of those results are presented in this section 
and the following one. Figure 8, for example, shows the streamlines and skin- 
friction variation, in physical coordinates, for a typical parabolic shear 

distribution. The relation between 
the physical and transformed 
c,7r / variables is 



C f 


2u /—— ■ X 
e / x 



Figure 8.- Streamlines for specified 
shear distribution. 


and 


Ip = 



( (50) 


dy 


/ 


In figure 9, the skin friction and 
streamline patterns for a different 
shear distribution are shown. For 
this case, the maximum reversed flow 
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Figure 9.- Streamlines for specified 
shear distribution. 
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Figure 10.- Velocity profiles for 
trailing-edge flow. 


occurs toward the reattachment side 
of the separation bubble. The divid- 
ing streamline has several rapid 
changes in slope, and this solution 
would be difficult to obtain if it 
were necessary to explicitly iterate 
for the location of the u = 0 line. 
Note that in all cases the normal 
coordinate is multiplied by the 
square root of the Reynolds number 
and that these solutions represent 
shallow separated regions confined 
to the interior of the viscous layer. 

The present method can also be 
used to calculate flows where 
reattachment occurs in a wake rather 
than on a solid surface. The details 
of this type of flow field in the 
immediate vicinity of the trailing 
edge are shown in figure 10. Here, 
the transition from boundary- layer 
flow to wake flow is assumed to occur 
on a scale that is small compared to 
the thickness of the viscous layer 
(see ref. 42). The prescribed 
boundary conditions of zero velocity 
and negative wall shear were thus 
discontinuously changed to zero 
shear and specified reversed-flow 
velocity at the trailing edge. Based 
on order-of -magnitude considerations, 
the initial reversed-flow velocity 
was taken to be equal to the value 
of the wall shear at the joining 


point. No attempt was made to ensure continuity of the dividing streamline 
or displacement thickness, although mass and momentum are conserved in the 
solution to the differential equations. 


Indications of Breakdown 

In the previous sections, it was demonstrated that the boundary- layer 
equations have regular solutions at separation and reattachment. The flow 
structure at the separation point agrees with the limiting form of the 
Navier-Stokes equations, and the Goldstein solution does not appear to be 
relevant for real flows. The square-root singularity in the boundary- layer 
equations is a consequence of specifying an external pressure distribution 
based on an inviscid solution determined as though there were no separation 
In practice, the pressure gradient is locally modified near the separation 
point such that the boundary- layer solution remains regular. The question 
that arises then concerns the manner in which the boundary- layer equations 
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eventually break down. Real flows tend to separate toward the rear of a closed 
body and vorticity is transported into the outer fluid. In some cases, the 
vorticity is confined to a relatively narrow region, or wake, downstream of 
the body. In other situations, behind a circular cylinder, for example, a 
large region of the fluid becomes rotational. The vorticity is no longer 
restricted to a thin viscous layer and the normal component of velocity 
ceases to be small compared with the' tangential component. In the present 
investigation, the region of separated flow is, of course, constrained to 
remain close to the surface, inside a layer of order 1//R. The numerical 
solutions may, however, suggest when this approximation is no longer realistic. 
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Figure 11.- Evidence of weak solutions 
for highly separated flow. 



Figure 12.- Normal velocity distribu- 
tion for highly separated flow. 


An indication of the possible 
breakdown in the boundary- layer 
equations is shown in figure 11 for 
a highly separated flow. As the 
mesh is refined, the computed values 
of m appear to become discontinuous 
at a point downstream of separation. 
Apparently, there are two solutions, 
one associated with separation and 
the other with reattachment, that 
are joined in the reversed-flow 
region. 

The distribution of v e , the 
transformed normal velocity, is 
shown on an expanded scale in fig- 
ure 12. The normal velocities 
increase rapidly downstream of the 
separation point, and the viscous 
layer begins to break away from the 
surface. Because of constraints 
imposed by the boundary conditions, 
however, a discontinuity in v (and 
in 3u/9x) occurs at the maximum 
value of v e , and the remaining 
solution is continuous. Although 
there is a certain degree of smooth- 
ing in the numerical results, the 
discontinuity in v e is evident in 
figure 12. A jump in v is an 
allowable weak solution of the dif- 
ferential equations and is apparently 
required for certain boundary condi- 
tions (e.g., large negative shears). 
If strong discontinuities occur when 
the shear distribution corresponding 
to a real flow is prescribed, how- 
ever, this can be taken to indicate 
the breakdown of the boundary- layer 
assumptions. 
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Figure 13.- Effect of shear variation 
in separated region. 
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Figure 14.- Effect of shear variation 
in vicinity of reattachment. 


The rapid variation of m, v e , 
and of the other flow quantities 
depends on the amount of reversed 
flow. This is illustrated in fig- 
ure 13 for a sequence of solutions 
where the specified shear distribu- 
tion was modified in the separated 
region. As the values of the shear 
become less negative, the solutions 
become increasingly smooth and con- 
tinuous. The streamlines correspond- 
ing to a = 0.1 were previously shown 
in figure 8. Even for this relatively 
mild case, the separating flow 
appears to undergo a rapid transition 
to the reattaching portion of the 
flow field at x = 2.7 approximately. 

The results of an additional 
numerical experiment are shown in 
figure 14. For this case, the wall 
shear was varied only in the down- 
stream portion of the separated 
zone and kept constant elsewhere. 

The nonlinearity and upstream influ- 
ence of the boundary- layer equations 
is evident in the computed distribu- 
tions of m and v e . Note also, how- 
ever, that the flow in the immediate 
vicinity of separation (x < 2.5) is 
not significantly affected by 
relatively large changes near 
reattachment. 


Upstream Influence 

Part of the success of 
approximate methods that use forward- 
marching schemes (e.g., refs. 9 and 
13) may be related to the limited 
upstream influence discussed above, 
particularly for flows with small 
separated zones. For the cases shown 
in figure 14, of course, it would 
not be possible to obtain accurate 
solutions downstream of x = 2.5 
without including the boundary con- 
ditions at reattachment. To investi- 
gate this question, calculations 
were made with the convective term 
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uu x set to zero for u ^ 0 with both the first- and second-order-accurate 
difference schemes used in a marching mode. Only backward differencing was 
employed for both, momentum and continuity, and the equations were completely 
relaxed at each x station before proceeding. 


The three point backward 
second-order scheme could be marched 
accurately a short distance into the 
separated zone. It always diverged 
rapidly, however, at approximately 
the location where m x became nega- 
tive. The first-order scheme, with 
moderate grid spacing, could be used 
for small bubbles but diverged for 
more separated flows. A typical 
calculation for a mildly separated 
flow is compared in figure 15 with 
an "exact" solution obtained using 
the correctly differenced second- 
order scheme with smaller step size. As the grid spacing was refined in x, 
the first-order marching began to diverge from the correct solution. The 
instability could be delayed by keeping Ay 5 Ax and by accepting a less 
. stringent iterative convergence criterion at each x station, but overall, 
the difference equations failed to converge to a solution as the grid was 
refined. 



Figure 15.- Comparison with forward- 
marching procedure. 


This experiment indicates that backward differencing, even with uu x = 0 
for u 5 0, is always unstable. For mild separation, the eigenvalues in the 
unstable range are small and dominate the numerical calculation only after a 
sufficient number of steps is taken. It is probable that the schemes of ref- 
erences 9 through 11 and 13 are also divergent, although they are useful for 
certain applications. 

To determine the effect of neglecting the upstream convection of momentum, 
additional calculations were performed with the term uu x set to zero for 
u < 0, but with the term u x in the continuity equation centrally differenced. 
In this manner, upstream influence is retained and the solution must again be 
obtained by relaxation methods. The results were essentially identical to the 
exact second-order solution, verifying that the upstream convection of 
momentum is not significant for laminar flows with limited separated regions. 


POSSIBLE EXTENSIONS 


An important extension of the present method is to match an inner, 
boundary- layer solution to an outer inviscid flow to calculate complete 
viscous-inviscid interactions. It would also be useful to compare results of 
the present method to experimental measurements of laminar separating and 
reattaching flows. Because low-speed boundary layers rarely remain laminar 
through reattachment, the computations must be extended to supersonic flows. 
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There are, for example, a number of reliable experiments for compression- 
corner interactions at supersonic speeds, as well as several different approxi- 
mate solutions and Navier-Stokes calculations available for comparison (e.g., 
refs. 43 and 44). It is indicated below how the method can be adapted to 
compressible flows, and an integral relation is proposed that offers promise 
of allowing the treatment of complete viscous-inviscid interactions. 


Compressible Flows 

To apply the method to compressible boundary layers, the following 
transformation can be used: 



If it is assumed that the density-viscosity product is constant through the 
layer, the following equations result: 


where 


- 3u _ 3u - ri — o*. ^ 3 2 u 
xu — + v — = m(l - u z ) + 


3x 


3y 


3y^ 


3u 3v u . , 

X 3^ + ^ + 2 11 + 

3y 


1 - yM 2 
.e 


1 + [(Y - D/2]M e 2 ; 


>m = 0 


x dM e 


u c 


m = — and M e = — 

M e dx e a e 


(52) 


These equations can then be solved in exactly the same fashion as equations 
(6a) and (6b), with M e calculated by integrating m. 


Viscous-inviscid Interations 

The solution for a complete interaction is complicated by the fact that 
x is specified. The following integral relation can, however, be used: 
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M 2 - 1 
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1 + [Cy - l)/2]M e 2 


/ 


u 2 dy 


m 


C53) 


For an assumed t distribution, the solution of equations (52) gives 
calculated values of m and hence of M e and p e . Using an inverse inviscid 
procedure, the distribution M e can be specified to obtain a new effective 
body shape, that is, the streamline slopes v e /u e . Then, from equation (53), 
a new estimate for x can be determined and the procedure continued until 
convergence is achieved. Based on recent experience with an integral scheme 
(ref. 8), it will probably not be advantageous to precisely match the 
intermediate iterations for v e /u e . 


It would, of course, be easier to specify v e directly for the viscous 
solution. For similar flows, an efficient scheme was developed by differen- 
tiating the continuity equation with respect to y and using standard second- 
order central differencing for v. The value of m was then updated by 
evaluating the continuity equation at the edge of the layer. This approach 
failed, however, for the complete boundary -layer equations with separated 
regions and was much slower for attached flows than the x specified schemes. 
An alternate approach, perhaps using the vorticity equation, may be required. 
All analytical and numerical evidence indicates, however, that the wall shear 
is the optimum boundary condition for calculating separated flows. 


CONCLUDING REMARKS 


The numerical procedures developed in the investigation provide an exact 
means for integrating the boundary- layer equations through separation and 
reattachment. The approach appears to be adaptable to the treatment of com- 
plete viscous-inviscid interactions for flow fields where the boundary layer 
remains confined to a narrow region: compression-corner flows or separation 

at the trailing edge of a streamlined body, for example. The method may also 
prove useful in evaluating different turbulence models for separated flows. 

As compared to complete Navier-Stokes solutions, the present approach allows 
an order-of-magnitude better resolution of the viscous region and requires 
considerably less computation time. Finally, a method based on the boundary- 
layer equations provides the most promising means for investigating the 
important problem of three-dimensional flow separation. 


Ames Research Center 

National Aeronautics and Space Administration 

Moffett Field, Calif., 94035, March 8, 1974 
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APPENDIX 


A program listing for the point-relaxation version of the first method 
is included. Only a description of the input and output is given; however, 
program variable names are the same as used in the text and should be self- 
explanatory. No effort was made to optimize the code or even to use a very 
efficient procedure for solving the attached region of the flow. A solution 
case corresponding to a = 0.1 in figure 13 is included. 


INPUT PARAMETERS 
(Subroutine INIT) 


JMAX = maximum number of points in x, 3 < JMAX < 120 
KMAX = maximum number of points in y, 4 - KMAX 5 100 
DX = Ax 

DY = Ay 

XO = x-location of initial profile 

UEXO = u e at XO 

SMO = m at XO 

ALPHU, ALPHV, ALPHM - relaxation parameters to update u, v, m 


(Subroutine PROFL) 

DYO = Ay in which initial profile is given 

KMAXO = number of data points to specify initial profile 

U(K),V(K) = u and v of initial profile 


(Main) 

ITERM = maximum number of iterations permitted 

RMAX = calculation is terminated if the maximum residual exceeds RMAX 

RMIN = residual at which iteration ceases and the converged solution is 
printed 

ALPHM2 = after an initial number of iterations, ALPHM is reset to this value 
ADDAL = increment to ALPHU and ALPHV after an initial number of iterations 

Wall Shear is analytically input in the present program. 


OUTPUT 


- The input parameters and the initial profile are printed. 

- Minimum output from a marching routine that calculates the attached flow 

region is printed. 
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Maximum residuals and their locations are printed every 10 iterations. 

The basic solution as a function of x is printed; data include j, x, m, 
u e , v e , and t. 

The solution profiles are printed at each x station. Data include k, y, 
u, v, y, u, v, and and interpolated values of y at constant values 

of ip/R. 



PROGRAM LISTING 
AND 

CASE RUN 
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PROGRAM FLOSEP ( INPUT , OUTPUT , TAPE5=INPUT , TAPE6=0UTPUT ) 


MAIM PROGRAM 

an ITERATIVE FINITE'DIFFERENCE METHOO FOR INTEGRATING the 
LAMINAR INCOMPRESSIBLE BOUNDARY-LAYER EQUATIONS THROUGH 
SEPARATION and reattachment 

COMMON SMC120) , SMC ( 1 20 ) » U ( 120 , 1 00 ) » V ( 1 20 , 1 00 ) , XB ( 1 20 ) 

COMMON /PARAM/ JMAX , KMAX , DX . DY , UEXO , 

1 ALPHU , ALPHV , ALPHM , RMAX , RM I N 

COMMON /RESID/ ITER , JREST , KREST , HESU , RESV,JU, KU 

C 

DIMENSION TAU ( 1 20 ) 

C 

C CONSTANTS 

ITERM = 1000 
RMAX = 10, 

RMIN s 0,00005 
ALPHM2 = 0,5 
ADDAL s 0 , 0 B 
C 

C INITIALIZATION 

CALL INIT 

SY2 = 1,/(DY*0Y) 

C 

C INPUT WALL SHEAR INTO TAU(J) ARRAY 

C EXAMPLE CASE 

ALF = 0.1 

TO = 0,33236/12, 

00 10 Jsl,JMAX 

X2 = X B ( J ) -2 , 

X6 s X2-R, 

IF ( X 2 ) 9,9,5 
5 IFCX6) 7,9,9 

7 TAUCJ) = T0*X2*X6*(1,+ALF*X2*X6) 

GO TO 10 

9 TAU(J) = T0*X2*Xfc 
10 CONTINUE 

C INITIALIZATION COMPLETE 

C 

C MARCHING IN ATTACHED FLOW REGIONS 
CALL MARCH CJ1, J2, TAU) 
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IF(J1'J2) 12,50,50 
12 CONTINUE 

relaxation part 

WRITE (6,500) 

C 

ITER = 0 
15 CONTINUE 
RMTST * 0,0 

JRM a 1 

C UPDATE M, EQ, 31 

DO 25 JaJl,J2 

RM a SM(J)+SY2*(O.<*U(J,2)-0,5»U(J,3)-3,*OY*TAU(J) ) 

SM( J) a SM(J)-4LPHM*RH 
SMC(J) a 0 , 5* ( SM ( J ) + l . 0 ) 

RMP a ABS (RM) 

IF(RMP-RMTST) 25,20,20 
20 RESM a RM 
RMTST a RMP 
JRM = J 
25 CONTINUE 

30 CONTINUE 

C CALL RELAXATION ROUTINE, METHOD ONE 

ITER » ITER +1 
CALL RELAX (Jl, J2 ) 

IK (ITER - ( ITER/ 1 0 ) * 1 0 ) 32,31,32 

C PRINT MAXIMUM residuals AND LOCATION 

C EVERY 10 ITERATIONS.. 

31 CONTINUE 

WRITE (6,501) ITER.RESV, JREST,KREST,RESU, JU,KU, RESM, JRM 

32 REST a ABS(RESV) 

IF (ITER 200) 38,30,38 

C CHANGE RELAXATION PARAMETER AFTER 200 ITER 

30 ALPHM S ALPHM2 

ALPHU a ALPHU + ADDAL 
ALPHV a ALPHV ♦ ADDAL 
IF (ALPHV - 1,0) 36,35,35 

35 CONTINUE 
ALPHU a 0.98 
ALPHV a 0,98 

36 CONTINUt 

WRITE(6,508) ALPHV, ALPHU, ALPHM 
C TEST WHETHER TO TERMINATE CALCULATION 

38 IF (REST-RMAX) 00,00,80 
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40 IF (ITER - ITERM) 45/45,42 
42 CONTINUE 

WHITE (6/502) 

GO TO 50 

45 IF (REST - RMIN) 46,46,15 

46 CONTINUE 
WRITE (6,507) 

C 

50 CONTINUE 

C TERMINATION WITH PRINT OUT 

CALL PRNT 
C 

60 CONTINUE 
STOP 

500 FORMAT (1H0,35X22HRELAX AT ION CALCULATION // 

1 7X»4HITER,5X,5HHES V,7X,7HJV KV,6X,5HRES U*7X, 

2 7H JU KU, 6X , 5HRES M,7X,2HJM) 

501 FORMAT (5X,I5,3(E15.5,215)) 

502 FORM AT ( 34H0 MAXIMUM ITERATIONS COMPUTED....) 

507 FORMAT (25H0 CONVERGED SOLUTION,,,.) 

508 FORMAT ( 1 HO , 1 7HALPH V , ALPHU , ALPHH ,3F13.5) 

END 

C 

SUBROUTINE init 

C INPUT DATA SUBROUTINE 

COMMON SM ( 1 20 ) , SMC ( 1 20 )., U ( 1 20 , 1 00 ) , V ( 1 20 , 1 0 0 ) , XB ( 1 20 ) 

COMMON /PaRAM/ JMAX , KMAX , DX , DY , UEXO , 

1 ALPHU , ALPHV , ALPHM , RMAX , RMIN 

DIMENSION Ul(lOO) , VI (100} 

C 

READ (5,500) JMAX , KMAX, DX, DY , XO , UEXO , SMO 
READ (5,501) ALPHU, ALPHV, ALPHM 

WRITE (6,505) JMAX, DX, XO, ALPHV, KMAX, DY, UEXO , ALPHU , SMO , ALPHM 

C INTERPOLATION OF INITIAL PROFILE IF NEEDED 

CALL PROFL (KMAX, DY, Ul, VI) 

DO 30 Jsl , JMAX 
XB(J) s XO+(J-l)*DX 
SM(J) = SMO 
SMC(J) * 0 , 5* (SMO+1 , ) 

DO 20 Kb 1 , KM AX 
U ( J , K ) b ui (K) 

V ( J, K ) s VKK) 

20 CONTINUE 
30 CONTINUE 
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RETURN 

500 FORMAT (2I5,5F10,0) 

501 FORMAT(BFIO.O) 

505 F0PMAT(1H1»35X»12H1NPUT VALUES// 

14X.6HJMAX a,lS,4X,4HDX a , F 1 0 . 5 , 6X , 4HX0 a , F 1 0 . 5 , «X , 7HALPMV a.FlO.S/ 
24X.6HKMAX 8, I5,4X,4HDY a,F 1 0 . 5 , 4X , 6HUEX0 a , FI 0 ,5 , «X , 7HALPNU a, 
3F10.5/ 39X.4HM0 a,F10.5,4X,7HALPHM «»F10.5) 

END 

C 

SUBROUTINE PROFL (KMAX1 , 0Y1, Ul, VI) 

C INTERPOLATION OF INITIAL PROFILE, 

DIMENSION Y(b0),U(60),V(60), C(5), S(5), T £5 ) , 

1 Yl(lOO), U1 ( 100) , VI Cl 00) 

INT c 2 

C INPUT INITIAL PROFILE 

READC5,501) OYO.KMAXO 
WRITE (6, 510) 

00 2 Kal.KMAXO 
READ(5,500) UCK),V(K) 

YIK) = (K-1)*DY0 

WRITE (6.511) X, Y ( K ) » U ( K ) » V(K) 

2 CONTINUE 

IF ( KMAX1 - KMAXO) 30,30,3 

3 CONTINUE 
KSAVE a 1 

DO 20 Kiel , KMA X 1 
Yl(Kl) a CK1-1 )*0V1 

4 DO 5 KaKSAVE, KMAXO 
KK a K 

IF (Y(K)-Yl (K1 ) ) 5,5,6 

5 CONTINUE 

6 IF(KSAVE-l) 9,9,7 

7 IFCY(KK-1)-Y1(K1)) 9,9,8 

8 KSAVE a 1 
GO TO 4 

9 KK = KK-tINT + n/2 
IF ( KK ) 10,10,11 

10 KK a 1 
GO TO 13 

11 M a KK+1NT 
IF(M-KMAXO) 13,13,12 

12 KK = KK-1 
GO TO U 

13 INTI = INT+1 
KSAVE a KK 

DO 1« L*1 , INTI 
C (L) = Y1 (Kl)-Y(KK) 

SCL) = U ( KK ) 

T(L) * V (KK ) 

14 KK a KK ♦ I 



DO 16 KK*1,INT 
I = KK+1 

15 D = C(KK)-C(I) 

SCI) s CC(KK)*S(I)-C(I)*SCKK))/D 
TCI) e CCCKK)*TCI)-CCI)*TCKK))/D 
I = HI 

IFCI-INT1) 15/15,16 

16 CONTINUE 
U1CM) = SCINT1) 

VICKI) s TCINT1) 

20 CONTINUE 
RETURN 

C NO interpolation 

30 DO 31 Ksl,KMAXl 
U1CK) a UCK) 

VI (K) = VCK) 

31 CONTINUE 

RETURN 

500 FORMAT (8F10.0) 

501 FORMATCF10,0,I5) 

510 FORM AT C 1 HO, 1 5 X,1«H INITIAL V ALUES//4X , 1HK , 7X , 1 HY , 1 1 X , 1 HU , 1 1 X , 1 HV ) 

511 F0RMATC2X,I3,3F12.6) 

end 

C 

SUBROUTINE MARCH CJ1, J2, TAU) 

C MARCHING IN ATTACHED REGION 

COMMON SMC120) , SMC(120),UC120,100),VC120,100),XBC120) 

COMMON /PARAM/ JMAX , KMAX , DX , DY , UEXO , 

1 ALPHU , ALPHV , ALPHM , RMAX , RMIN 

DIMENSION TAUC120) , UX C 1 00 ) 

TAUWT b 0.02 
ALPHM2 s ,5 
SY a . 5*DY 
SYY = DY*DY 
SX2 a ,5/DX 
SY2 * l./CDY*DY) 

C 

N 3 2 
JINT » -l 
KM a KMAX-1 

IF CTAUC3) - TAUNT) 50,50,5 
C 

5 N s N +1 
J = N 

J1 a J t JINT 
J2 a J 
JINT = 0 
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IF ( J * JMAX) 6,6,50 


C TEST TO SEE IF PROFILE IS ATTACHED 

6 IFCTAUCJ) - TAUNT ) 50,50,7 

7 CONTINUE 

C OBTAIN GOOD GUESS BY USING EXTRAPOLATION OF LAST COMPUTED PROFILES 
DO 6 K=2,KM 

U ( J, K) b 2,*U(J-1,K) -U(J-2,K) 

8 V ( J , K ) « 2,*V(J-1,K) - V(J-2,K) 

SM ( J ) a SM(J-l) 

c 

ITER * 0 
UXC1) = 0. 

C 

10 ITER = ITER + l 
DO 20 Jaj 1 , J2 
JR s J-l 
JRR = J-2 
REST » 0. 

RM a SMCJ)+SY2*(«,*U(J,2)-0.S*U(J,3)-3.*DY»TAU(J)) 

SM C J 3 = SM(J)-ALPHM*RM 
SMC ( J) = 0.5*CSM(J)+t,0) 

DO 18 K = 2 , KM 
KR b K-l 
KP a K *1 
C 

IF ( J-2) 12,12,13 

12 UX ( K ) a 3X2* t U(J + 1,K) - UCJR,K)) 

U2X a UX(K) *( UCJ+l.K) ♦ U(JR.K)) 

DIAX a 0. 

GO TO 14 

13 UX(K) » SX2* ( 3,*U(J,K) -<t.*U(JR,K) +U(JRR,K)) 

U2X a SX2*( 3.*U(J,K)**2 -a. *U < JR, K ) **2 ♦ U(JRR,K)**2) 

DIAX = 3 , *XB ( J) *SX2 
1« CONTINUE 

UY a SY*(U(J,KP)-UCJ,KR) ) 

FU = UY*V(J/K)*SYY*C0.5*XBCJ)*U2X-SM(J)*(1.-U(J>K)**2)) 

RU a U(J,KR)-2.*UCJf K)+U(J,KP)-FU 

RV a VCJ,K)-VlJ,KR)+SY*(XB(J)*CUXCK)fUX(KR))*SMCCj)*CUCJ*K)* 

1 U(J,KR))3 

OT a 2. + SYY*U(J,K)*CIAX 
OU a RU/DT 
DV a - HV 

U ( J ,K ) a U ( J , K ) + OU 
V C J , K ) a V ( J , K 3 ♦ DV 
RT a ABS(RV) 

IFCR7 - REST) 18,18,15 
15 REST = RT 
18 CONTINUE 
20 CONTINUE 
C 

IF C ITER - 20) 26,26,25 
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25 ALPHM s ALPHM2 

26 CONTINUE 

IF < REST -RMAX) 27,100,100 

27 IF ( REST - RMIN) 30,30,28 

28 I F ( ITER - 600) 10,100,100 

30 CONTINUE 

IF ( N"3 ) 35 , 35, RO 

35 WRITE(6,501) 

12 = 0 
RZ = 0, 

00 36 J = 1,2 

VtJ,KM A X) s V(J,KM) - OY*SMCCJ) 

TAuw = .5*1 -3,*u(0»l) ♦R.*U(J,2) -U(J»3))/DY 
WRITE (6, 500) J, IZ, RZ,VCJ, KMAX), TAUW, SMCJ) 

36 CONTINUE 
«0 J a N 

TAUW * , 5* ( -3 , *U ( J, 1 ) YR.oUCJ.Z) -UCJ,3))/DY 
VCJ.KMAX) * V(J,KM) - DY*SMC(J) 

WRITE (6, 500) J, ITER, REST, VCJ,KMAX), TAUW, SM(J) 

GO TO 5 
C 

50 J1 = N 
J2 s UMAX 
RETURN 

100 J s N 

TAUW = ,5*( “3, *U C J, 1 ) +«.*UCJ,2) -UCJ,3))/0Y 
V(J,KMAX) = V(J,KM) - DYaSMCCJ) 

WR I TEC 6, 50 0) J, ITER, REST, VCJ,KMAX), TAUW, SM(J) 

STOP 

500 FORMAT C 1H ,2I5,4F13.5) 

501 F0RMAT(1H1, 35X18HMARCHING PROCEDURE // 5X1MJ, 2X4HITER, 

1 7X5HRES V, 8X4HVMAX, 1 0 X3HT AU , 1 0 X l HM ) 

END 

C 

SUBROUTINE RELAX CJ1, J2) 

C RELAXATION FOR REVERSED FLOW 

COMMON SMC120) , SMC ( 1 20 ) , U C 1 2 0 , 1 0 0 ) , V ( 1 20 , 1 0 0 ) , XB ( 1 20 ) 

COMMON /PaRAM/ JMAX , KMAX , DX OY , UEXO , 

• I ALPHU , ALPMV , ALPHM , RMAX , RMIN 

COMMON /RESID/ ITER , JREST , KREST , RESU , RESV,JU, KU 

DIMENSION UXC100) 

c 

EPS1 = .015 
EPS2 * 0.005 
KM s KMAX -1 
C 

EPS = EPS1-EPS2 
RESU = 0. 
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RESV a 0.0 
RTEST a o, 

REST * 0.0 
SY a 0 , 5*DY 
SYY = DY*DY 
SX2 a 0.5/DX 
00 SO JnJl . J2 
JR a J-l 
JP a J+l 
UX(1) a 0.0 
DO 30 Ka 2, KM 
KR a K.l 
KP = Kfl 
C 

c 

IF(J-JMAX) 10,13,13 

10 UX ( K ) a SX2«(UIJP,K)-UCJR,K)) 

IFCJ-2J 11,11,12 

11 U2X = SX2»(U(JP,K)**2-U(JR,K)**2) 

D 1 AX « 0. 

GO TO 20 

12 T a UCJ,K) 

IFCT-EP31) 16,16,10 

13 UX ( K ) a SX2*(3,*UCJ«K>-4,*UCJR»K)+UCJR-1,K)) 

C 

C ATTACHED FLOW 

C 

14 JQ a J.2 ' 

U2X a SX2*(UCJQ,K)**2-4,*UCJR,K)**2+3.*U(J,K)*«2) 

D I A X = 3.*XB(J)*SX2 
GO TO 20 
C 

C SEPARATED FLOW 

c 

16 U2X a SX2*(U(JP,K)**2-U(JR,K)**2) 

D 1 AX a 0. 

IF (T*EPS2) 20,20,18 
C 

C SEPARATION POINT 

C REATTACHMENT POINT 

C 

18 JO a J-2 
U2P a U2X 

U2X » SX2*(U(JQ,K)**2-4,*U(JR,K)**2+3,*U( J,K)**2) 

TA a (T-EPS2J/EPS 
U2X * TA*U2X+(1,-TA)»U2P 
DIAX a 3,*TA*XB(J)*SX2 
C 

20 CONTINUE 

UY a SY*(U(J,KP)-UCJ,KR)) 

FU a UY*V(J,K)+SYY*(0.S*XB(J)*U2X-SM(J)»(1.-U(J,K)*«2)) 

RU a U(J,KRW,*UCJ,K)+U(P,KP)-FU 

RV a VCj,K)-VlJ,KR)tSY*(X8(J)*(UX(K)tUX(KR))FSMCtJ)*CU(J»K)t 
1 U(J,KR») 

DT = 2. ♦ SYY*U(J,K)»DIAX 





DU a NU/DT 
DV : • BV 

UCJ,K) a U(J,K)+ALPHU*DU 
VCJ,K) a V(J,K)+ALPHV*DV 

RT a ABS(RU) 

IF (RT»RTEST ) 22,22,21 

21 RESU a RU 
RTEST s RT 
JU a J 

KU a K 

22 CONTINUE 

RTT = ABSCRV) 

IF (RTT-REST ) 27,27,26 

26 RESV a RY 
REST a RTT 
JREST a J 
KRE S T a K 

IF (REST-RMAX) 27,27,100 

27 CONTINUE 
30 CONTINUE 

K a KMAX 

RV = V(J,KMAX)-V(J,KM)+SY*(XBCJ)*UXCKM)+SMC(J)*(1.YU(J»KM))) 
V(J,KMAX) a V(J,KMAX)-ALPHV*RV 
RTT = ABSCRV) 

IF (RTT-REST ) 35,35, 34 
3a RESV * RV 
REST * RTT 
JREST a J 
KREST a K 
35 CONTINUE 
50 CONTINUE 
100 RETURN 
END 

SUBROUTINE PRNT 


C OUTPUT SUBROUTINE 


COMMON SMC120) 

COMMON /PARAM/ JMAX 
1 ALPHU 

COMMON /RESID/ iter 
COMMON /STRM/ PS I (100) 
1 KIM 


, SMC(120),U(120,100),V(120,100),XB(120) 

, KMAX , DX , DY , UEXO , 

, ALPHV , ALPHM , KMAX , RM I N 

, JREST , KREST , RESU , RESV,JU, KU 
, YSTCIOO) , POUT (60 ) , YK60) , y2(60) , 
, K2M 


DIMENSION UY ( 120) , UE(120) , F(120) , FX ( 1 20 ) , USTC100) , 
1 VST(IOO), ETA(IOO) 

HY a 0.5/DY 
KEND a KMAX 
SX a 0.5*DX 
DO 10 J=1 ,JMAX 

UY(J) = HY*(-3.*UCJ,l)+a,*U(J,2)-UCJ,3)) 

10 CONTINUE 
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REST * RESV 

WRITE 16,500) ITER, JR£ST,KREST, REST 
C 

C INTEGRATION FOR UECX) 

C 

UE ( 1 ) a UEXO 
F(l) = ALOG (UEXO ) 

IF C XB ( 1 ) -.00001 ) 12,13,13 

12 XB(1) a ,00001 

13 CONTINUE 

F X ( 1 ) * SMCD/XBCl) 

00 15 Ja2, JMAX 
FX(J) a SM(J)/XB(J) 

F ( J) a F(J-1HSX»(FX(J)+FX(J-1)) 

UE(J) a EXPCF(J)) 

15 CONTINUE 

WRITE (6,502) 

WRITE (6,503) (J,XB(J)»SM(J),U£(J)«V(J*KMaX)»UY(J)»Js1,JMAX) 

DO 20 Kal , KMAX 
ETA(K) a (K-l ) *DY 
20 CONTINUE 
YST(l) * 0, 

VST ( 1 ) a 0. 

PSI(l) a 0. 

DO 50 Ja 1 , JMAX 

Cl a SQRT(XB(J)/UE(J)) 

C2 a l./Cl 

C3 a SQRT(XB(J)*UE(J)) 

51 a 0.5»(SM(J)-1.) 

52 a 0 ,5*DY*C3 
UST(l) a UE(J)*U(J,1) 

DO 30 Ka2» KMAX 
Y a ETA(K) 

YST(K) a C 1 *Y 
UST(K) a U(J,K)*UE(J) 

VST(K) a (V(J,K)-S1»Y*U(J,K) )*C2 
PSI(K) = PSI(K-1)+S2*(U(J,K)6U(J,K-D) 

30 CONTINUE 

DSTR a YST(KMAX)-PSI(KMAX)/UE(J) 

WRITE(6,510) XB(J) .DSTR 
C 

CALL STREAM (KENO) 

C 

WRITE (6,511) 

IF (KIM ,EQ. 0) GO TO «0 

WRITE (6,512) (K,ETA(K),U(J,K),V(J,K),YST(K),UST(K), VST(K), 

1 PSl(K), POUT (K) , Yt(K), Y2(K), Kal, KIM) 

C KIM - NUMBER OF POINTS WITH SEPARATION 


«0 CONTINUE 



K2 a KIM + 1 

WRITE (6/514) (K,ETA(K),U(J,K),V(J,K),YST(K),UST(k), VST(K), 

1 PSICK), POUT (K)/ Y 1 C K ) » KSK2/K2M) 

C K2M - TOTAL NUMBER OF INTERPOLATED POINTS, 

C (SEPARATED PLUS SINGLE VALUED) . 


K3 a K2M ♦ 1 

WRITE (6,513) (K,ETA(K),U(J,K),V(J,K),YST(K),UST(K), VST(K), 

1 PSI(K) , KaK3» KMAX ) 

50 CONTINUE 
RETURN 

500 FORMAT ( 9h 1 ITER a , 1 5 , 9H JRESV a,I5,9H KRESV s,I5,8h RESV a, 
1 M3. 5) 

502 FORMAT(5HO J,6X,4MX(J),eX,2HSM,9X,6HU EDGE / 7X, 6HV EDGE / 6X , 

1 10HDU/DY WALL) 

503 FORMAT (2X,I3,F10,3,4F13.6) 

510 F0RMAT(1H1,6HX(J) a F 1 2 .5 , 2X 8HDELSTR = F12.6 ) 

511 FORMAT (t 02 X, 12H INTERPOLATED/ 
13X,1HK,4X,3HETA,9X,1HU,10X,1HV,15X,1HY,11X,3HUST,9X,3HVST,9X, 

2 3HPSI,15X,3HPSI,7X,1HY/9X,1HY) • 

512 FORMAT(I4,F8.3,2F12,6,4X,4F12.6,IOX,F7.4,2F10.5) 

513 FORMAT (I4,F8.3,2F12.6,4X,4F12.6) 

514 FORMAT (I4,F8.3,2F12.6,4X,UF12.6,10X,F7.4,FlO t 5) 


C 


END 

SUBROUTINE STREAM (KMAX) 


C INTERPOLATION FOR STREAM FUNCTION 


COMMON /STRM/ PSKlOO) / YST(IOO) / P0UT(60) , Yl(60) » Y2(60) » 
1 KIM , K2M 

DIMENSION C ( 4 ) , S ( 4 ) 


PSMIN a -0.10 
PSMAX a 10,0 

NO a »200,*PSMIN 
NO a NOW 1 
N 1 a NO+9 
N2 a Nl+10 
N3 a N2+10 
DO 8 Ns 1 , 60 
IF (N-N3) 2,1,1 

1 POUT ( N ) o 6+N-N3 
NMAX a N 

IF (POUT(N)-PSMAX) 8,9,9 

2 IF(N-N2) 4,3,3 

3 POUT(N) a 0.5*(N-N2+1) 
GO TO 8 
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9 IFCN-N1) 6,6,5 

5 POUT (N) a O,0b*(N-Nl) 

GO TO S 

6 POUT (N) s 0,005* ( N-NO ) 

8 CONTINUE 

9 CONTINUE 

c 

C**»**FIND MINIMUM PSI 
INT = 2 

DO 10 Ka2,KMAX 
K K - K 

IF(PSl(KK)-PSI(KK-m 10,10,11 

10 CONTINUE 

11 KPMIN a KK-1 
IF(KPMIN-INT) 12,12,15 

12 NMIN s NO 
KPMIN a 1 
GO TO <10 

15 PM I N s PS I (KPMIN) 

C 

C*****FIND INITIAL PRINTOUT VALUE 

DO 17 N? 1 , NMAX 
NN = N 

IF(POUTCN)-PMIN) 17,17,18 

17 CONTINUE 

18 NMIN = NN 

N 1 a NO-NMIN 
C 

C*****INTERPOLATION.,.PSI from wall TO U a 0 

DO 30 L c 1 , N 1 
NN a NO-L 
DO 20 K=l, KPMIN 
KK a K 

IF(POUTINN)-PSICKK)) 20,20,21 

20 CONTINUE 

21 KK a KK-ClNT + n/2 
IF(KK) 22,22,23 

22 KK a 1 
GO TO 25 

23 M a KK+INT 
IF(M-KPMIN) 25,25,29 

29 KK a KK-1 
GO TO 23 

25 INTI a INT+1 
DO 26 J=1 , INTI 

C t J ) = POUTCNN)-PSI(KK) 

SCJ) = YST(KK) 

26 KK a KK ♦ 1 

DO 28 J a l,INT 

I = J + l 

27 SCI) = (C(J)*SCI)-C(I)*5(J)J/(C(JJ-C(I)) 
I a Itl 



IF(I-INTt) 27,27,28 
26 CONTINUE 

Y2(NN) * S(INTl) 

30 CONTINUE 

KPMIN a KPMIN+t 
00 CONTINUE 

KSAVE a KPMIN 
C 

C»****INTERPOL*TION, , ,PSI FROM U = 0 TO EDGE 

DO 60 NeNMIN.NMAX 
NN a N 

DO 05 k«ksa ve , kma x 

K K — K 

IF(PS1(KK)-P0UT(NN)) 05,05,06 
05 CONTINUE 

4,4 IF (KK-KPMIN) 07,07,08 

07 KK = KPMIN 
GO TO 09 

08 M a KK+2 
IF(M-KMAX) 09,09,06 

09 KSAVE a KK 
DO 50 jet, 3 

C(J) a POUT(NN)-PSI (KK) 

S(J) a VST (KK ) 

50 KK a KK+1 

S(2) = (C(1)*S(2)-C(2)*S(1))/(C(1)-C(2)) 

S (3) a (C(1)*$(3)-C(3)*S(1))/(C(1)-C(3)) 

T l (NN) a (C(2)*S(3)-C(3)*S(2))/(C(2)-C(3)) 

60 CONTINUE 
K1 a 0 

IF(NO-NMIN) 75,75,65 
65 NOt = NO-1 

DO 70 L=NMIN,NOl 

Kt a K1 ♦ t 

POUT (K 1 ) a POUT (L) 

Y 1 ( K 1 ) a Tl (L) 

Y2(Kt) a Y2(L) 

70 CONTINUE 
75 CONTINUE 
KIM a Kl 
K2 = KIM 
DO 80 L a N0,NMAX 
K2 a K2 + 1 
P0UT(K2) 8 POUT (L) 

Y t (K2) a Y 1 ( L ) 

80 CONTINUE 
K2M a M2 
RETURN 

END 

52 52,125 .20 0, t, 0, 
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25 50 

0,00000 0 , 

0,08304 - 0 . 

0,16597 - 0 . 

0.24848 . 0 . 

0.33003 > 0 . 

0,40991 - 0 , 

0.48726 - 0 . 

0.56111 » 0 . 

0,63047 > 0 , 

0.69442 - 0 . 

0.75216 - 0 . 

0,60313 - 0 . 

0,84704 . 0 . 

0,88390 . 0 , 

0,91400 > 0 , 

0,93790 - 1 . 

0.95632 - 1 . 

0,97010 - 1 , 

0,98010 - 1 . 

0,98712 - 1 . 

0,99191 - 1 , 

0,99506 - 1 . 

0,99708 - 1 , 

0,99832 - 2 . 

0,99906 - 2 , 

0,99949 - 2 , 

0.99973 - 2 . 

0,99986 - 2 . 

0,99993 - 2 . 

0,99996 - 2 , 

0,99998 - 2 . 

0,99999 - 3 . 

0,99999 - 3 , 

0,99999 - 3 , 

0,99999 - 3 . 

0,99999 - 3 . 

0,99999 - 3 . 

0,99999 . 3 , 

0,99999 . 3 , 

0,99999 . 4 , 

0,99999 - 4 , 

0.99999 - 4 , 

1,00000 > 4 . 

1,00000 . 4 , 

1,00000 - 4 , 

1,00000 - 4 . 

1,00000 . 4 , 

1,00000 . 5 . 

1,00000 - 5 , 

1,00000 . 5 . 

JENDDS 


00000 

00519 

02075 

04665 

08281 

12906 

18513 

25065 

32513 

40793 

49834 

59555 

69869 

80687 

91924 

03498 

15337 

27377 

39566 

51861 

64230 

76649 

89100 

01571 

14055 

26546 

39041 

51538 

64037 

76537 

89036 

01536 

14036 

26536 

39036 

51536 

64036 

76536 

89036 

01536 

14036 

26536 

39036 

51535 

64035 

76535 

89035 

01535 

14035 

26535 



INPUT VALUES 


JMAX =52 DX s .12500 

UMAX a 52 DY s .20000 



INITI 

AL VALUES 


K 

Y 

U. 


1 

0.000000 

0,000000 


2 

.250000 

.083040 


3 

.500000 

.165970 


4 

.750000 

,248480 


5 

1.000000 

,330030 


6 

1.250000 

.409910 


7 

1,500000 

,487260 


6 

1,750000 

.561110 


9 

2.000000 

,630470 


10 

2.250000 

,694420 


1 1 

2.500000 

.752160 


12 

2.750000 

.803130 


13 

3.000000 

.847040 


14 

3,250000 

.883900 


15 

3.500000 

.914000 


16 

3.750000 

.937900 


17 

4.000000 

.956320 


18 

4.250000 

,970100 


19 

4,500000 

.980100 


20 

4.750000 

.987120 


21 

5.000000 

.991910 


22 

5.250000 

.995060 


23 

5.500000 

,997080 


24 

5.7S0000 

.998320 


25 

6.000000 

.999060 


26 

6.250000 

.999490 


27 

6.500000 

.999730 


28 

6.750000 

,999860 


29 

7.000000 

.999930 


30 

7.250000 

,999960 


31 

7,500000 

,999980 


32 

7,750000 

,999990 


33 

8,000000 

,999990 


34 

8.250000 

.999990 


35 

8.500000 

.999990 


36 

8,750000 

.999990 


37 

9,000000 

,999990 


38 

9 , 2i> 00 0 0 

.999990 


39 

9,500000 

.999990 


40 

9.750000 

.999990 


41 

10.000000 

,999990 


42 

10.250000 

,999990 


43 

10,500000 

1,000000 


44 

10.750000 

1,000000 


45 

11.000000 

1.000000 


46 

11.250000 

1.000000 


47 

11,500000 

1,000000 


48 

1 1.750000 

1.000000 


49 

12.000000 

1,000000 


50 

12.250000 

1,000000 



X0 = 0,00000 ALPHV * 

UEXO = 1.00000 ALPMU = 

MO = 0.00000 ALPHM 3 


V 

0.000000 
.005190 
.020750 
.046650 
.082810 
, 129060 
. 185130 
.250650 
.325130 
.407930 
.498340 
.595550 
.698690 
.806870 
.919240 
1.034980 
1.153370 
1.273770 
1,395660 
1.518610 
1.642300 
1,766490 
1.891000 
2.015710 
2. 140550 
2.265460 
2.390410 
2,515380 
2.640370 
2.765370 

2.890360 

3.015360 

3.140360 

3.265360 

3.390360 
3.515360 
3.640360 
3.765360 

3.890360 

4.015360 

4.140360 

4.265360 

4.390360 
4,515350 
4,640350 
4,765350 
4.890350 
5,015350 
5.140350 
5,265350 


,90000 

,90000 

,05000 
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MARCHING PROCEDURE 


iKarntr 3 X 3 XXXXX 33 


3 3 KIMK 1 f\J 


ruru— *«-»©3r^4!4>4>4>4>4)xxxxxx 
pnmpnpomrurururururururupurufurururu 


OBSo^a O 3 <h® Mir>- 3 rn ro 
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X(J) ■ 2.62500 DEISTR e 7.525652 

INTERPOLATED 

K ETA U V V - UST VST PSJ PSI Y Y 

1 0,000 0,000000 0,000000 0,000000 0,000000 0,000000 0.000000 -.0150 2,42487 2,12359 
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X £ J 3 = 2,75000 DELSTR s 8.043320 

INTERPOLATED 

K ETA U V Y UST VST PSI PSI Y Y 

1 0 ■ 000 0,000000 0.000000 0,000000 0.000000 0,000000 0,000000 -.0200 3,08916 1,64859 
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X(J) ■ 2,87500 DEL8TR » 8,<<87083 
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